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Abstract. Simulations of past climates require altered 
boundary conditions to account for known shifts in the Earth 
system. For the Last Glacial Maximum (LGM) and subse- 
quent deglaciation, the existence of large Northern Hemi- 
sphere ice sheets caused profound changes in surface to- 
pography and albedo. While ice-sheet extent is fairly well 
known, numerous conflicting reconstructions of ice-sheet to- 
pography suggest that precision in this boundary condition is 
lacking. Here we use a high-resolution and oxygen-isotope- 
enabled fully coupled global circulation model (GCM) 
(GISS ModelE2-R), along with two different reconstructions 
of the Laurentide Ice Sheet (LIS) that provide maximum and 
minimum estimates of LIS elevation, to assess the range of 
climate variability in response to uncertainty in this bound- 
ary condition. We present this comparison at two equilibrium 
time slices: the LGM, when differences in ice-sheet topog- 
raphy are maximized, and 14ka, when differences in maxi- 
mum ice-sheet height are smaller but still exist. Overall, we 
find significant differences in the climate response to LIS 
topography, with the larger LIS resulting in enhanced At- 
lantic Meridional Overturning Circulation and warmer sur- 
face air temperatures, particularly over northeastern Asia and 
the North Pacific. These up- and downstream elfects are as- 
sociated with differences in the development of planetary 
waves in the upper atmosphere, with the larger LIS resulting 
in a weaker trough over northeastern Asia that leads to the 
warmer temperatures and decreased albedo from snow and 
sea-ice cover. Differences between the 14ka simulations are 
similar in spatial extent but smaller in magnitude, suggesting 


that climate is responding primarily to the larger difference 
in maximum LIS elevation in the LGM simulations. These 
results suggest that such uncertainty in ice-sheet boundary 
conditions alone may significantly impact the results of pa- 
leoclimate simulations and their ability to successfully sim- 
ulate past climates, with implications for estimating climate 
sensitivity to greenhouse gas forcing utilizing past climate 
states. 


1 Introduction 

The Last Glacial Maximum (LGM, ~ 21 ka) provides a valu- 
able target to test the ability of general circulation models 

(GCMs) to simulate a climate for which they were not de- 

signed (Mix et al., 2001; 1PCC, 2007; Braconnot et al., 2012). 
In addition, the large changes in temperature and greenhouse 
gases from the LGM to present may provide a possible means 

of assessing climate sensitivity to changes in atmospheric 
CO 2 (Crucifix, 2006; Hansen et al., 2008; Schmittner et al., 
2011, 2012; Fyke and Eby, 2012; Hargreaves et al., 2012; 
PALAEOSENS, 2012). The last deglaciation (~20 to 7ka) 
was the most recent period in which changes in Earth’s or- 
bit around the Sun caused Northern Hemisphere ice-sheet re- 
treat (Clark et al., 2009; Carlson and Winsor, 2012; He et al., 
2013) and rising atmospheric greenhouse gas concentration 
(Monnin et al., 2001; Lemieux-Dudon et al., 2010), which 
provides an evolving climate state (Shakun and Carlson, 
2010) for testing GCMs (e.g., Timm and Timmerman, 2007; 
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Liu et al., 2009; Shakun et al., 2012; Gregoire et al., 2012; 
He et al., 2013). For instance, by 14 ka, Northern Hemi- 
sphere ice sheets were still relatively large (60-70 % remain- 
ing by area; Dyke, 2004; Gyllencreutz et al., 2007), but at- 
mospheric greenhouse gas concentrations had already risen 
by ~60% of their total deglacial change (Monnin et al., 
2001; Lemieux-Dudon et al., 2010) and boreal summer inso- 
lation was ~ 7.5 % higher than present/LGM levels (Berger 
andLoutre, 1991). 

However, the simulated glacial state in any one model is 
sensitive to the boundary conditions used as a starting point 
for the simulation (Manabe and Broccoli, 1985; Broccoli and 
Manabe, 1987; Hyde and Peltier, 1993; Justino et al., 2006; 
Abe-Ouchi et al., 2007; Liu et al., 2009; Pausata et al., 2011; 
Hofer et al., 2012). Thus, uncertainty in a particular set of 
glacial boundary conditions may overshadow a GCM’s sim- 
ulated change in climate and its ability to constrain climate 
sensitivity using the LGM as a starting point. Along this line, 
Hargreaves et al. (2012) suggested that disregarding high- 
latitude response at the LGM avoids bringing glacial bound- 
ary condition ambiguity into climate sensitivity work, al- 
though this omission then ignores many climate feedbacks. 
Some of the important boundary conditions are well docu- 
mented, such as the orbital parameters that drive the magni- 
tude and seasonality of insolation (Berger and Loutre, 1991), 
as well as the concentration of atmospheric greenhouse gases 
(Monnin et al., 2004; Jouzel et al., 2007; Ltithi et al., 2008; 
Lemieux-Dudon et al., 2010). Other boundary conditions of 
past glacial climates have proven to be more elusive. The 
aerosol forcing (particularly dust and black carbon) at the 
LGM and through the deglaciation is known only at in- 
dividual locations, with limited spatial resolution (Petit et 
al., 1981; Thompson et al., 1995; Mahowald et al., 1999, 
2006a, b; Reader et al., 1999; Power et al., 2008; Ganopolski 
et al., 2010). While large uncertainties still remain on the 
impact of aerosols on radiative forcing (Penner et al., 2004; 
Lohmann and Feichter, 2005; Forster et al., 2007; Chylek and 
Lohmann, 2008), the LGM dust loading may issue a forcing 
on the order of ~ — 1 W m -2 that is comparable to changes 
in greenhouse gases alone (-■ — 2.85 W m -2 ), at least in the 
tropics, where ice-sheet albedo is not a factor (Harrison et al., 
2001; Claquin et al., 2003; Crucifix, 2006). Likewise, recon- 
structions of LGM and deglacial vegetation and the related 
impact on surface albedo are limited to low-resolution global 
approximations and higher resolution estimates that are only 
regional in scope (Jackson et al., 2000; Prentice and Jolly, 
2000; Harrison et al., 2001; Ray and Adams, 2001; Bigelow 
et al., 2003; Williams, 2003; Pickett et al., 2004; Williams 
et al., 2011). The impact of such uncertainties in boundary 
conditions may even be amplified through possible internal 
feedbacks of climate models with the inclusion of new dy- 
namic components, such as vegetation and dust (Ridgwell 
and Watson, 2002; Mahowald et al., 2006a, b; O’ishi and 
Abe-Ouchi, 2013). 


Perhaps the greatest source of known uncertainty in glacial 
boundary conditions relates to ice-sheet thickness. Although 
the geographical extents of LGM and deglacial ice-sheets are 
fairly well mapped (Denton and Hughes, 1981, 2002; Dyke 
and Prest, 1987; Anderson et al., 2002; Clark and Mix, 2002; 
Bennike and Bjorck, 2002; Dyke, 2004; Gyllencreutz et al., 
2007), direct observations of ice thickness and topography 
are limited and usually absent for the highest/thickest por- 
tions of the ice sheets (e.g., Denton and Hughes, 1981, 2002; 
Dyke et al., 2002; Clark, 1992; Clark et al., 2009; Goehring 
et al., 2008; Carlson and Clark, 2012). Therefore, ice-sheet 
height must be simulated through geophysical or glaciolog- 
ical modeling approaches. The reconstructions of ICE-4G, 
ICE-5G, and ICE-6G (Peltier, 1994, 2004; Toscano et al„ 
2011) have proven to be useful ice-sheet boundary condi- 
tions to employ in GCMs due to their geophysically based 
solutions, global scope, and product availability. However, 
the differences between these reconstructions (Peltier, 2004; 
Toscano et al., 2011), as well as other ice-sheet-specific 
reconstructions (Clark et al., 1996; Licciardi et al., 1998; 
Tarasov and Peltier, 2004; Tamisiea et al., 2007; Argus and 
Peltier, 2010; Lambeck et al., 2010) suggest that there is a 
large range in existing ice-sheet boundary conditions, partic- 
ularly for the Laurentide Ice Sheet (LIS). 

This large uncertainty between reconstructions of ice- 
sheet geometry is unfortunate considering that glacial to- 
pography is one of the dominant drivers of atmospheric and 
oceanic circulation, temperature, and precipitation in simula- 
tions of deglacial climate (Justino et al., 2006; Abe-Ouchi et 
al., 2007; Pausata et al., 2011; Hofer et al., 2012; Tharammal 
et al., 2013). Glacial orography leads to reduced surface air 
temperatures through vertical lapse rate alone, but the pres- 
ence of these large ice mountains also has been shown to have 
downstream effects, altering upper-atmosphere flow (Cook 
and Held, 1988; Bromwich et al., 2004; Justino et al., 2006; 
Abe-Ouchi et al., 2007; Langen and Vinther, 2009), midlat- 
itude storm tracks and wintertime precipitation (Kageyama 
et al., 1999; Hofer et al., 2012), and Atlantic Meridional 
Overturning Circulation (AMOC) strength (Yu et al., 1996; 
Adkins et al., 2002; Wunsch, 2003; Curry and Oppo, 2005; 
Lynch-Steiglitz et al., 2007; Otto-Bliesner et al., 2007; Weber 
et al., 2007; Arzel et al., 2008). 

Here we test the impact of LIS geometry on LGM and 
deglacial climate using two alternate reconstructions of LIS 
topography that provide upper and lower bounds for this 
boundary condition. As an upper bound, we use ICE-5G 
(VM2) (Peltier, 2004; hereafter 5G), which has the high- 
est LIS topography of any LIS reconstruction, with a dom- 
inant ice dome over the western Keewatin sector (Fig. 1). 
As a lower bound, we employ the dynamics-driven flow- 
line reconstruction of Licciardi et al. (1998, hereafter L) 
(Fig. 1). The maximum height of this reconstruction is com- 
parable to ICE-4G (Peltier, 1994; Clark et al., 1996), but it 
is ~ 20 % lower than that of ICE-5G, particularly over the 
western LIS (Fig. 1). The use of these reconstructions is 
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Licciardi et al. Topography, 21 ka 
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Fig. 1. Laurentide Ice Sheet topographies used in 21 ka (top row panels) and 14 ka (bottom row panels) simulations. The left column presents 
the Licciardi et al. (1998) reconstructions used in 21ka-L and 14ka-L. The middle column presents the ICE-5G reconstructions (Peltier, 2004) 
used in 21ka-5G and 14ka-5G. For comparison, the single plot in the right column shows the LIS topography used in the PMIP3 boundary 
conditions (http://pmip3.lsce.ipsl.fr/) 


a particularly relevant assessment of the range of variabil- 
ity in the Paleoclimate Modelling Intercomparison Project 3 
(PMIP3) topographic boundary conditions, as the maximum 
ice-sheet heights of the three ice-sheet reconstructions used 
to create the averaged topography of the LIS in the PMIP3 
(http://pmip3.lsce.ipsl.fr/) lie in between the upper and lower 
bounds used in our study. We compare the effects of these 
LIS boundary conditions at two time periods, 21 and 14 ka, 
to assess the impact of varying LIS topography under dif- 
ferent climate forcing conditions. The 21 ka simulations pro- 
vide a test of topographic variability under full LGM bound- 
ary conditions with different LIS maximum elevations. At 
14 ka, both LIS topographies are reduced and differences be- 
tween maximum LIS elevation are smaller (~15%). The 
14ka simulations provide a test of ice-sheet geometry dif- 
ferences under a different forcing scenario, with changes in 
precession, obliquity, and increased atmospheric greenhouse 
gas concentrations. 

2 Methods 

2.1 Model 

The NASA Goddard Institute for Space Studies (GISS) Mod- 

elE2, coupled to the Russell (R; Russell et al., 2000) ocean 

is a fully coupled atmosphere-ocean global climate model 
(Schmidt et al., 2014). We use the same version of ModelE2- 

R as is being used for the Coupled Model Intercomparison 

Project Phase 5 (CMIP5) simulations (NINT physics ver- 
sion), with an atmosphere resolution of 2° latitude by 2.5° 


longitude with 40 vertical layers up to 0.1 mb, and an ocean 
resolution of 1° latitude by 1.25° longitude with 32 depth 
layers. The 21 ka simulations presented here are the GISS 
ModelE2-R contribution to LGM experiment of CMIP5 (en- 
sembles r 1 i lp 1 50 and rli lp 151; see Sect. 2.3). ModelE2- 
R includes passive water isotopologue tracers (i.e., H) 2 * * * * * 8 0, 
HDO, hereafter referred to as water isotopes), removing 
one degree of uncertainty when comparing its diagnostics 
to paleoclimate reconstructions based on water isotope vari- 
ability (e.g., ice cores, speleothems, ocean foraminifera) to 
test model skill (e.g., LeGrande et al., 2006; Carlson et al., 
2008a, b; LeGrande and Schmidt, 2009). 

2.2 Boundary conditions 

Each simulation uses appropriate insolation for the time pe- 
riod owing to changes in orbital parameters (Berger and 
Loutre, 1991), along with appropriate glacial greenhouse gas 
concentrations (Joos and Spahni, 2008) (Table 1). Coastline 
and basin geometry were adjusted to reflect an ~ 120 m low- 
ering in sea level at 21 ka and a ~ 86 m lowering at 14 ka, 
in accordance with the volume of water contained in ice- 
sheet boundary conditions and glacial isostatic adjustment 
(Peltier, 2004). The volume of water held in the two alternate 
LIS reconstructions is different, but in order to avoid small 
coastal idiosyncrasies between the two simulations for each 
time slice, we kept the ocean bathymetry (and the land-sea 
mask) the same for each pairing of LIS boundary conditions. 
Because we aim to study differences in ocean circulation due 
to changes in the LIS topography alone, bathymetry was held 
consistent between the simulations at each time slice so as 
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Table 1 . Boundary conditions for experiments presented in this paper. The abbreviations for the orbital parameters refer to eccentricity (E), 
obliquity in degree (O), and angular precession (longitude of perihelion, or omega) in degree (P) (Berger and Loutre, 1991). 


Simulation 

Time slice 

Orbital 

parameters 

Greenhouse gas 
concentrations 

Sea 

level 

change 

LIS geometry and 
max elevation 

21 ka-L 

21 ka 

E: 0.019398 
O: 22.989 
P: 113.98 

CO 2 : 188 ppm 
CH 4 : 385 ppm 
N 2 O : 200 ppm 

-120 m 

Licciardi et al. (1998) 
3560 m 

21ka-5G 

21 ka 

E: 0.019398 
O: 22.989 
P: 113.98 

CO 2 : 188 ppm 
CH 4 : 385 ppm 
N 2 O : 200 ppm 

-120 m 

Peltier (2004) 
4520 m 

14ka-L 

14 ka 

E: 0.020180 
O: 24.004 
P: 228.37 

CO 2 : 239 ppm 
CH 4 : 630 ppm 
N 2 O: 261 ppm 

— 86 m 

Licciardi et al. (1998) 
2890 m 

14ka-5G 

14 ka 

E: 0.020180 
O: 24.004 
P: 228.37 

CO 2 : 239 ppm 
CH 4 : 630 ppm 
N 2 O: 261 ppm 

— 86 m 

Peltier (2004) 
3400 m 

Control 

Pre industrial 
(Oka) 

E: 0.017236 
0 : 23.446 
P: 101.37 

CO 2 : 285 ppm 
CH 4 : 791 ppm 
N 2 O: 275 ppm 

0 m 

Modem 


not to introduce basin differences along ocean boundaries 
that might impact ocean dynamics. In addition, ModelE2-R 
ocean uses straits - essentially pipelines for ocean tracers and 
mass - from present day to accommodate sub-grid-scale pas- 
sages. All straits were eliminated (because of sea level low- 
ering) except for Gibraltar and Bab Al Mandab (Red Sea), 
which were reduced, and no new straits were created. 

ModelE2 does not have a dynamic vegetation component, 
so terrestrial vegetation cover was assigned according to 
mapped conditions at the LGM (25 000-15 000 yr BP; Ray 
and Adams, 2001) for regions not covered by ice sheets. 
A separate reconstruction for 14 ka does not exist; there- 
fore we continued to use this LGM reconstruction, modified 
slightly to account for changes in ice-sheet extent by expand- 
ing nearby vegetation types into these areas no longer cov- 
ered by ice. Continental river routing was assigned in accor- 
dance with ice-sheet configuration and its impact on drainage 
basins (e.g., Licciardi et al., 1998, 1999; Peltier, 2004). 

2.3 Ice-sheet geometry sensitivity 

To test the influence of ice-sheet topography on glacial cli- 
mate, two separate simulations were conducted at each time 
slice using two alternate reconstructions of LIS topogra- 
phy: one using the geophysically based ICE-5G (5G; Peltier, 
2004; CMIP5 ensemble rlilpl50 for the LGM), and the 
other using the same extent of ICE-5G but replacing the LIS 
sector topography with that of an alternate reconstruction 
based on a flow-line model that simulates glacier dynamics 
over deformable and rigid beds with varying till viscosities 
(L; Licciardi et al., 1998; CMIP5 ensemble rlilpl 5 1 for the 


LGM). We use the maximum reconstructions from Licciardi 
et al. (1998) that incorporate higher specified effective till 
viscosities. These Licciardi et al. (1998) maximum recon- 
structions were concluded to provide a better representation 
of the LIS than the “minimum” reconstructions due to better 
agreement with other ice-sheet reconstructions and relative 
sea level records (Clark et al., 1996; Licciardi et al., 1998). In 
both simulations, all other ice-sheet topographies and extents 
(Greenland, Scandinavian, etc.) are assigned using ICE-5G. 
We focus on the LIS because it was the largest of the Qua- 
ternary ice sheets with presumably the greatest effect on cli- 
mate (Alley and Clark, 1999; Clark et al., 1999; Clark and 
Mix, 2002; Carlson and Clark, 2012). 

The two LIS reconstructions have significantly different 
topographies (Fig. 1). At 21 ka, the 5G reconstruction pro- 
vides a much taller LIS, particularly the Keewatin Dome 
over western Canada, with elevations across this region at or 
above 4000 m. The 21 ka-L reconstruction is lower in max- 
imum height and moves more of the ice mass to the east 
with three distinct ice domes centered over Keewatin, south- 
ern Ontario, and central Quebec. Comparatively, the 5G-LIS 
has a lower topography over eastern Ontario and Quebec, but 
the large Keewatin Dome dominates the western topography. 
Additionally, the 5G reconstruction has more abrupt changes 
in topography, whereas the L reconstruction has smoother 
transitions from high to low elevations, as would be expected 
from a flow-line model of ice deformation over hard and 
soft beds and similar to inferences from the geologic record 
(Dyke and Prest, 1987, Jenson et al., 1996). 
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At 14 ka, both reconstructions retain the same general 
geometry as 21 ka, but absolute topography is diminished 
(Fig. 1). ICE-5G still has a dominant Keewatin dome, while 
the L reconstruction continues to have three prominent 
domes. The L reconstruction exhibits only a small diminu- 
tion in ice topography over the Ontario and Quebec domes, 
such that a west-east dipole in topography difference exists 
between the two reconstructions (Fig. 1). 

Our climate model simulations were started before the 
PM1P3 protocol for LGM conditions had been determined, 
but our selection of LIS boundary conditions provide lim- 
its of topographic uncertainty that go into the PMIP3 LIS 
elevations. The LIS in the PMIP3 boundary conditions is 
similar in geometry to ICE-5G (with higher Keewatin dome 
over western Canada and lower topography over eastern On- 
tario and Quebec), but the maximum topography of the LIS 
in PMIP3 is more similar to the LIS reconstruction used in 
21ka-L (Fig. 1). 

2.4 Equilibrium simulations - 21, 14, and 0 ka 

A control simulation (Oka) with water isotopes was initiated 
from a 700 yr long spinup from Levitus and Boyer (1994a, b) 
and Levitus et al. (1994) and run an additional 1000 yr. The 
21 and 14ka simulations were also initiated from the same 
700 yr long spinup from Levitus and Boyer (1994a, b) and 
Levitus et al. (1994), but these experiments have altered 
mean ocean <5 18 0, <5D, and salinity (see section below), the 
last of which requires a much longer integration for the 
ocean density structure to adjust. Surface air temperature 
(SAT) drift is significant for the first 500 yr in each of the 
deglaciaFglacial simulations - some of which is related to a 
temporary enhancement of AMOC due to the method with 
which enhanced salinity was applied. Fractional increases in 
salinity were uniformly added to each grid box, where in 
reality additional salinity would have likely been preferen- 
tially buried in the deep ocean. Such colder and saltier con- 
ditions are more conducive to the enhanced overturning seen 
in the initial spinup, which has been shown to require an ex- 
tended initialization (Zhang et al., 2013). After 1500 yr of 
integration, SAT drift is less than 0.03 °C per century, sur- 
face ocean temperature drift less than 0.02 °C per century, 
and deep ocean temperature drift less than 0.03 °C per cen- 
tury and in each of the simulations, but overall drift is never 
absent. Note that these simulations run 500 yr longer than 
those present in the CMIP5 archive to insure that the water 
isotopologues have come into equilibrium (water isotopes are 
not a part of the CMIP5 archive). Initially, each simulation 
had an enhanced ocean circulation, which diminished with 
approach toward equilibrium. In addition, cold surface tem- 
peratures drove sea ice growth to the bottom of the ocean 
grid in some shallow coastal polar regions, requiring an ex- 
pansion of the land mask as such “ice shelves” are defined 
in other regions of the boundary conditions. Here we present 
the final 100 yr of these simulations. 


2.5 Water isotopes 

ModelE2-R includes water isotopologue tracers that allow 
for the direct comparison of <5 18 0 measurements taken from 
ice cores, speleothems, ocean sediments, and other pale- 
oclimate records (Carlson et al., 2008a, b; LeGrande and 
Schmidt, 2009). We have selected terrestrial data that pro- 
vide a proxy of the /> 18 0 in precipitation (hereafter <5 18 O a ) 
throughout the LGM and deglacial climate (see Table SI in 
the Supplement for list of <5 18 0 proxy records used in this 
analysis). To estimate LGM and deglacial <5 18 0 in the ocean 
(hereafter <$ 18 0 0 ), data from ocean sediment cores with con- 
tinuous coverage from the LGM to the preindustrial era were 
selected (Table SI in the Supplement). We used the average 
of all data within 1000 yr of each time slice (i.e., 20-22 ka 
for 21 ka, 13-15 ka for 14 ka, 2-Oka for control) to estimate 
the (5 18 O a and <5 18 0 0 anomaly from each record for direct 
comparison with the model. 

3 Results 

Our 21 and 14 ka simulations are broadly consistent with pre- 
vious GCM simulations of LGM and deglacial climate (e.g., 
CLIMAP, 1981; Kageyama et al., 1999; Justino et al., 2004; 
Otto-Bliesner et al., 2006, 2007; Abe-Ouchi et al., 2007; 
Braconnot et al., 2007a, b; Timm and Timmerman, 2007; Liu 
et al., 2009). However, this paper is primarily concerned with 
the differences that arise due to variation in the LIS topogra- 
phy alone. Since most of the climate differences directly over 
the ice sheet are due to orographic differences alone, we fo- 
cus on the LIS impacts in other regions in the Northern Hemi- 
sphere where GCM anomalies are maximized. For a detailed 
description of modeled 21 and 14ka anomalies relative to the 
G1SS ModelE2-R preindustrial control simulation, please see 
the Supplement. 

3.1 Atmosphere differences 

3.1.1 Surface air temperature 

SAT anomalies are significantly different between the two 
simulations at 21 ka (Fig. 2a). Global area-weighted mean 
SAT anomalies (21 ka minus control) for the 21ka-5G and 
21ka-L simulations are —5.1 ± 0.1 °C and —5.4 ± 0.1 °C, re- 
spectively (unless otherwise stated, uncertainty is expressed 
as standard deviation across decadal variability about the 
100 yr mean). This difference in SAT is significant, devel- 
oping solely in response to changes in LIS boundary con- 
ditions. Since differences in SAT between our two scenar- 
ios are exacerbated directly over the LIS due to topographic 
differences (i.e., vertical lapse rate of temperature and an 
~ 1000 m maximum LIS height difference), the global mean 
SAT anomalies with the LIS area removed are —4.7 ± 0.1 °C 
and — 5.0 ± 0. 1 °C, thus confirming that differences between 
the simulations are not due to differences directly over the 
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Jet Speed Difference, 21 ka 500mb height Difference, 21 ka 



Fig. 2. Differences between 21 ka simulations (21ka-L minus 21ka-5G) for the following annually averaged atmospheric variables: (a) surface 
air temperature (°C); (b) precipitation minus evaporation (mm day -1 ); (c) <5 18 O a (%o); (d) surface wind speed (ms -1 ); (e) atmospheric jet 
speed (defined as the wind speed at 250 hPa, ms -1 ); (f) geopotential height at the 500 mb level, with zonal mean removed (m); (g) ground 
albedo (%); and (h) planetary albedo (%).The LIS is masked out to focus on downstream differences. Ice-sheet extents outlined in bold black 
lines. 


ice sheet alone, but rather due to differences in other regions 
that significantly impact the global mean (Fig. 2a). 

The most prominent downstream differences in SAT be- 
tween the 21 ka simulations occur over northeastern Asia, 
Beringia, and the North Pacific, where 21ka-L is 6-9 °C 
colder than 21ka-5G (Fig. 2a). Proxy records from this region 
are limited, but one recent modern analog approach based on 
pollen spectra from Lake ETgygytgyn in northeastern Siberia 
suggests that the mean temperature of the warmest month 


(July) at 20 ka is ~6°C colder than present (Melles et al., 
2012), whereas the model results in a cooling of July tem- 
peratures by 16.5 ± 1.1 °C (21ka-L) and 12.8 ± 0.7 °C (21ka- 
5G). Model results appear to be much too cold in this region, 
but the lake record may not extend far enough back to cap- 
ture full LGM cooling, and lake ice permanence at the LGM 
may have strongly limited pollen transport into the lake dur- 
ing this time, biasing the record towards warm years (Melles 
et al., 2007, 2012). One SAT reconstruction from fossil lipids 
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in paleosols of eastern Asia suggests LGM summer cooling 
of 8-10°C (Peterse et al., 2011), whereas our 21 ka sim- 
ulations suggest summer cooling of 9.0±0.7°C (21ka-L) 
and 8.1 ± 0.4 °C (21ka-5G). Reconstructions from Zagoskin 
and Burial lakes in northwestern Alaska suggest July LGM 
cooling of 4 ± 1.6 °C and 5 ± 1.6 °C, respectively, which is 
consistent with general LGM cooling of the region (Viau et 
al., 2008; Kurek et al., 2009). Our 21 ka simulations suggest 
July LGM cooling of 8.8 ± 0.9 °C (21ka-L) and 4.6 ± 0.4 °C 
(21ka-5G) for Zagoskin Lake, and 12.2 ± 1.0 °C (21ka-L) 
and 8.5±0.7°C (21ka-5G) for Burial Lake. Other pollen- 
based reconstructions may actually provide some evidence 
for slightly warmer in Alaska during the LGM (Bartlein et 
al., 2011), although there may be no-analog issues in this re- 
gion due to the limited number of individual records avail- 
able in this synthesis. Unfortunately, no sites from northeast- 
ern Asia exist in the LGM climate reconstructions of Bartlein 
et al. (2011), limiting the model-data comparison of SAT, 
where our 21 ka simulations provide the greatest differences. 

For the 14ka simulations, SAT anomalies are similar 
for the two different LIS topographies, with global mean 
SAT anomalies of —3.0 ± 0.1 °C (14ka-L) and —2.8 ± 0.1 °C 
(14ka-5G). However, regional distinctions still exist between 
the simulations. Cooling over northeastern Asia and the 
North Pacific continues to occur in the 14ka-L simulation rel- 
ative to 14ka-5G by 1-2 °C (Fig. 3a). Model results suggest 
a July SAT change in this location of — 0.8±0.7°C (14ka- 
L) and ±0.6±0.9°C (14ka-5G), consistent with the record 
at Lake El’gygytgyn suggesting July temperatures near mod- 
ern conditions at 14 ka (Melles et al., 2012). The other proxy 
records from eastern Asia and northwestern Alaska suggest 
SAT near present conditions (±1.6 °C; Peterse et al., 2011; 
Viau et al., 2008; Kurek et al., 2009). Our model results 
indicate colder 14 ka anomalies of — 2.7±0.5°C (14ka-L) 
and -2.0 ±0.8 °C (14ka-5G) for eastern Asia, -2.0 ± 0.6 °C 
(14ka-L) and — 1.5±0.7°C (14ka-5G) for Zagoskin Lake, 
and —2.6 ± 0.7 °C (14ka-L) and -1.9±0.7°C (14ka-5G) 
for Burial Lake. These small negative SAT anomalies may 
suggest that the model simulations are slightly too cold in 
this region relative to proxies that suggest anomalies closer to 
zero (±1.6°C; Peterse et al., 2011; Viau et al., 2008; Kurek 
et al., 2009), but the model results are still consistent with the 
proxies within uncertainty. Like with the LGM simulations, 
however, the lack of proxy records across this region makes 
such model-data comparisons tenuous (Bartlein et al., 201 1). 

3.1.2 Precipitation minus evaporation 

The location of the Intertropical Convergence Zone (ITCZ) 
is consistent between 21 ka-L and 21ka-5G, but 21ka-L simu- 
lates greater precipitation north of the Equator in the tropical 
Pacific (Fig. 2b). The 21ka-5G simulation shows greater pre- 
cipitation over the northern mid-latitude Pacific and across 
eastern America, south of the LIS, indicative of a southward 


displacement of the primary storm tracks (Justino et al., 
2004). 

At 14 ka, differences in P — E between simulations are 
small, but highlight a slight enhancement of Pacific precipi- 
tation south of the Equator in 14ka-5G (Fig. 3b). Small dif- 
ferences in precipitation also occur over the North Atlantic, 
with drier conditions in 14ka-L along the east coast of North 
America and the Atlantic coast of Europe and slightly wetter 
conditions in the interior of the North Atlantic. 

3.1.3 6 ls O of the atmosphere differences 

Differences in /> 18 O a largely reflect the differences in SAT, 
with the greatest relative depletion of 21ka-L occurring in 
northeastern Asia and the North Pacific, where SATs show 
the greatest cooling relative to 21ka-5G (Fig. 2c). In the 14 ka 
simulations, there is a decoupling between the differences in 
SAT and the differences in <5 18 O a , with relatively enriched 
<5 18 O a in 21ka-L across northeastern Africa and into south- 
central Asia (Fig. 3c). Additionally, the 21ka-L shows rela- 
tively enriched values across much of the Arctic, but depleted 
values immediately downwind of the LIS and over southern 
Greenland. 

3.1.4 Atmospheric circulation 

The Northern Hemisphere subpolar jet in 2 lka-5G is stronger 
and more zonal across the mid-latitude Pacific and to the 
south of the LIS (Fig. 2e). This difference in jet location is 
associated with the increased precipitation in the 21ka-5G 
across eastern North America as the primary storm tracks are 
confined further to the south than in 21ka-L (see Sect. 3.1.2; 
Fig. 2b). The subpolar jet in 21ka-L is more wavelike and 
lies to the north of the 21ka-5G jet across most of Asia and 
the Pacific. Differences in surface wind speed reflect this jet 
displacement, particularly over the North Pacific, where the 
more northern 21ka-L jet results in greater surface winds. 
Across much of the North Atlantic, 21ka-5G surface winds 
are stronger than 21 ka-L, where both simulations are already 
enhanced relative to Oka (Fig. Sid in the Supplement). 

Differences in 500 mb height between the two 21 ka sim- 
ulations express a shift in the location and depth of the 
dominant stationary wave patterns, particularly over Siberia 
and Beringia, where 21 ka-L heights are substantially lower 
than those in the 21ka-5G simulation (Fig. 2f). These lower 
heights in 21ka-L reflect a deepening of the trough over east- 
ern Asia and a weakening of the ridge over Beringia relative 
to 21ka-5G. This general reduction of 500 mb heights across 
this region provide a greater influence of Arctic air masses 
that help to drive the colder SAT (Sect. 3.1.1). 

The higher LIS in 14ka-5G still induces an increase in jet 
speed and southward displacement of the polar and subtrop- 
ical jets relative to 14ka-L (Fig. 3e). However the stronger 
subpolar jet in 14ka-5G now extends further to the north over 
northern Europe. The general differences in 500 mb heights 


www.clim-past. net/1 0/487/201 4/ 


Clim. Past, 10, 487-507, 2014 


494 


D. J. Ullman et al. : Assessing the impact of Laurentide Ice Sheet topography on glacial climate 


SAT Difference, 14 ka P-E Difference, 14 ka 



A 18 Oa Difference, 14 ka Surface Wind Speed Difference, 14 ka 




Fig. 3. Same as Fig. 2 but with the atmospheric differences between 14 ka simulations (14ka-L minus 14ka-5G). 


and the implied stationary wave pattern in the 21 ka simula- 
tions still exists between the 14 ka simulations, with a contin- 
ued lowering of 14ka-L 500 mb heights in northeastern Asia 
and Beringia that is associated with the colder SAT in this 
region relative to 14ka-5G (Fig. 3f). 

3.1.5 Albedo 

The 21ka-L simulation presents a region of higher albedo in 
northeastern Asia, where there is additional snow cover over 
Siberia and Beringia, along with a relative expansion of sea 
ice over the Sea of Okhotsk (Fig. 2g). The 21ka-5G has a 
higher albedo in a region along the southern margin of the 


LIS, likely due to the southward tracking of the subpolar 
jet and enhanced snowfall. The disparity in the area extent 
of these albedo changes results in globally averaged ground 
albedo of 18.5 % for 21ka-L compared with 18.1 % for 21ka- 
5G. Planetary albedo differences mimic the ground albedo 
differences (Fig. 2h), although cloud cover over the region 
in northeastern Asia is enhanced by 8-10% in the 21ka-L 
simulation (not shown). 

At 14 ka, surface albedo is largely the same between 
the two simulations, except for the expansion of sea ice 
over the Nordic Seas in 14ka-5G, which drives an increase 
in albedo relative to 14ka-L (Fig. 3g). Globally averaged 
ground albedo is 16.3 % for 14ka-L and 16.4 % for 14ka-5G. 
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SST Difference, 21 ka 
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Atlantic Stream Function difference, 21 ka 
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Fig. 4. Differences between 21 ka simulations (21ka-L minus 21ka-5G) for the following annually averaged ocean variables: (a) sea surface 
temperature (°C), (b) sea surface salinity (psu), (c) <5 1S 0 0 (d) sea ice fraction (%), (e) Atlantic Ocean overturning stream function (Sv), and 
(f) Pacific Ocean overturning stream function (Sv). 


Again, planetary albedo differences mimic ground albedo 
differences, suggesting the prominence of surface albedo 
changes in overall albedo (Fig. 3h), with cloud cover differ- 
ences of less than 2 %. 

3.2 Ocean differences 

3.2.1 Ocean temperature differences 

Averaged globally, SST anomalies are slightly colder 
in 21ka-L (— 2.7±0.1°C) compared with 21ka-5G 
(— 2.5 ± 0.1 °C). This 21ka difference extends throughout 
the entire ocean with total global mean ocean temperature 
anomalies of-2.15 ± 0.01 °C (21ka-L) and -2.01 ± 0.01 °C 
(21ka-5G) (global mean ocean temperatures are scaled by 
mass of ocean water). Greater global cooling in 21ka-L 
persists in the deep ocean, where temperature anomalies 
averaged over the bottom 2000 m are — 1.18 ±0.01 °C 
(21ka-L) and -1.04 ± 0.01 °C (21ka-5G). The global mean 


ocean anomaly for the 21ka-L simulation overlaps with 
Kr/N 2 estimates of global mean ocean cooling at the LGM 
of 2.7 ± 0.6 °C (Headly and Severinghaus, 2007). 

Regional differences in SST between the 21ka-L and 
21ka-5G simulations correlate with SAT differences (Figs. 2a 
and 4a), particularly in the North Pacific, Sea of Okhotsk, 
and Bering Sea, where 21ka-L simulates colder SST by 
up to 3.5 °C. The large differences between the simulations 
in the North Pacific indicate that this is a sensitive region 
for testing the model results. SST proxy records are sparse 
from the North Pacific, but limited data seem to indicate 
colder conditions. Alkenone and transfer function SST re- 
constructions from off the coast of Oregon document LGM 
cooling of 4-7 °C (Doose et al., 1997; Ortiz et al., 1997; 
Mangelsdorf et al., 2000; Flerbert et al., 2001; Rosell-Mele 
et al., 2004), which are more consistent with the 21ka-L re- 
sult of 4.0±0.3°C cooling than with the 21ka-5G cooling 
of 3.1 ± 0.2 °C. Alkenone records from the Sea of Okhotsk. 
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Fig. 5. Same as Fig. 4 but with the ocean differences between 14 ka simulations (14ka-L minus 14ka-5G). 


however, suggest LGM SST cooling of only 0-1 °C (Harada 
et al., 2006; Seki et al., 2004, 2007). The 21ka-5G simulation 
has LGM cooling of 4.6 ± 0.2 °C in this region, and the 21ka- 
L simulation shows cooling of 7.4 ± 0.2 °C, suggesting that 
neither simulation adequately captures SST in this somewhat 
isolated basin. 

At 14 ka, globally averaged SST anomalies are equiva- 
lent between the two ice-sheet topographies (—1.5 ± 0.1 °C 
in 21ka-L and — 1.4±0.1°C in 21ka-5G), but total global 
mean ocean anomalies are still colder in the 14ka-L sim- 
ulation with an anomaly of — 0.83 ±0.1 °C (21ka-L) rela- 
tive to — 0.71 ±0.1 °C (21ka-5G). In the deep ocean (bot- 
tom 2000 m), however, temperature is nearly consistent 
with Oka with anomalies of — 0.09 ±0.01 °C (14ka-L) and 
±0.05 ± 0.01 °C (14ka-5G). Differences between the 14 ka 
simulations are most pronounced across the North Pacific 
(Fig. 5a), where 14ka-L SSTs are colder, similar to SAT dif- 
ferences over these regions (Fig. 4a). 


3.2.2 Sea surface salinity differences 

Surface waters are generally fresher in 21ka-L, particu- 
larly in the Arctic Ocean and the Gulf of Mexico, reflect- 
ing a greater contribution of LIS melt water through the 
MacKenzie and Mississippi rivers (Fig. 4b). The 21ka-5G 
simulation presents a greater contribution of Scandinavian 
Ice Sheet melt water through the Ob and Yenesei rivers, lead- 
ing to a localized freshening of sea surface salinity (SSS) 
along the coast, but in general, Arctic surface waters are more 
saline in 21ka-5G. Northern tropical Pacific surface waters 
are slightly fresher in 21ka-L in association with enhanced 
P — E over this region (Fig. 2b). 

SSS differences at 14 ka are largely similar to those in the 
21 ka simulations, with fresher waters in 14ka-L in the Arctic 
Ocean, the Gulf of Mexico, and along the Atlantic coast of 
North America (Fig. 5b). The 14ka-L simulation continues to 
present localized freshening due to greater melt water contri- 
butions of the LIS to the MacKenzie and Mississippi rivers. 
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SSS differences are negligible over the Pacific, reflecting the 
minimal changes in P — E (Fig. 3b). 

3.2.3 5 ls O of the surface ocean differences 

The differences in <5 18 0 0 at 21 ka are largely similar to the 
differences in SSS, with depleted values in 21ka-L relative 
to 21ka-5G across the Arctic Ocean and the Gulf of Mexico 
(Fig. 4c), reflecting the enhanced melt water contributions in 
21ka-L to the MacKenzie and Mississippi rivers (Fig. 4b). In 
addition, 21ka-L waters are depleted along the north shore 
of the Scandinavian Ice Sheet, consistent with reduced SSS 
and enhanced melt water contribution in 21ka-5G along this 
margin. The Sea of Japan shows slightly depleted <5 I8 0 0 val- 
ues in 21ka-L but no change in SSS, suggesting a shift in 
the isotopic composition of runoff arriving to this basin but a 
consistent volume of this runoff between the simulations. 

In the 14ka simulations, however, this relationship be- 
tween SSS and <5 18 0 0 becomes decoupled, with fresher arctic 
conditions in 14ka-L associated with more enriched values 
of <5 18 0 0 relative to 14ka-5G (Fig. 5c). However, the 14ka- 
L simulation continues to show depleted values in relation 
to the greater contribution of fresh water from the MacKen- 
zie and Mississippi drainage outlets. In addition, the Sea of 
Japan continues to reflect a difference in the isotopic compo- 
sition of river runoff to the basin, but here the 14ka-5G sim- 
ulation presents depleted <5 18 0 0 values. Again, the lack of a 
SSS difference in this basin suggests that freshwater runoff 
is consistent between the simulations (Fig. 5b). 

3.2.4 Sea ice differences 

Annually averaged sea ice differences are negligible across 
much of the ocean basins, except in the Sea of Okhotsk and 
Bering Sea, where the 21ka-L has a significant increase in sea 
ice over that of 21ka-5G (Fig. 4d). The expansion of sea ice 
across the Nordic Seas in the 21 ka simulations relative to the 
control simulation is consistent between 21ka-L and 21ka- 
5G, with 21ka-L having a slightly greater sea-ice fraction (no 
more than 5 %). 

In the 14 ka simulation, global sea ice differences are small 
(Fig. 5d). Differences in sea ice are minimized in the Sea of 
Okhotsk compared to the 21 ka simulation. The only notable 
differences occur in the Nordic Seas, with enhanced sea ice 
in the 14ka-5G simulation. 

3.2.5 Ocean circulation differences 

All simulations result in a mean AMOC that is stronger 
than the preindustrial control simulation (see Supplement; 
Oka control AMOC strength: 28.2 ± 0.7 Sv). The 21ka-5G 
simulation exhibits a mean AMOC transport that is ~ 10% 
stronger than 21ka-L (33.2 ± 0.7 Sv and 30.8 ± 0.6 Sv, re- 
spectively). The difference in AMOC stream function reflects 
not only the increase in AMOC transport for 21ka-5G but 
also a deepening of the main overturning and a slight increase 


in Antarctic bottom water in flow (Fig. 4e). In the Pacific, 
the 21ka-5G simulation shows an enhanced contribution of 
Antarctic intermediate and bottom waters throughout most 
of the basin (Fig. 4f), with the 21ka-5G having the great- 
est Antarctic bottom water influence (i.e., deepwater stream 
function more negative than 21ka-L). In addition, the 21ka- 
5G shows enhanced production of North Pacific intermedi- 
ate water relative to 21ka-L (or less of a reduction relative to 
Oka). 

Kinematic and water-mass proxy records continue to re- 
fine reconstructions of AMOC at the LGM, with overturn- 
ing strength anywhere from 30 to 40% weaker than present 
(McManus et al., 2004) to about the same as today (Lynch- 
Stieglitz et al., 2007; Ritz et al., 2013), and even the pos- 
sibility of more rapid overturning at the LGM (Gherardi 
et al., 2009; Lippold et al., 2012). Water mass tracers sug- 
gest the shoaling of glacial North Atlantic deepwater and a 
greater contribution of Antarctic bottom water (Curry and 
Oppo, 2005), but this increased stratification may also im- 
ply a more vigorous AMOC. Such uncertainty in circula- 
tion strength and depth is mirrored in a diverse array of 
modeled AMOC results, some of which present stronger 
AMOC during the LGM (Otto-Bliesner et al., 2007), simi- 
lar to our results. In the Pacific, the expansion of enriched 
<5 13 C throughout the deep ocean suggests the increased influ- 
ence of Antarctic-sourced waters in this basin (Matsumoto et 
al., 2002; Herguera et al., 2010). Other North Atlantic simu- 
lations have shown that enhanced AMOC is associated with 
the strengthening of the deep ocean circulation in the Pacific 
(Chikamoto et al., 2012), which is consistent with the en- 
hanced negative stream function in our 21 ka results (Fig. S3f 
in the Supplement). 

The 14ka simulations also present strong AMOC rela- 
tive to the preindustrial control (see Supplement). The 14ka- 
5G simulation continues to have stronger AMOC transport 
(32.7 ±0.7 Sv for 14ka-5G; 30.5 ± 0.6 Sv for 14ka-L), but 
the depth of overturning is similar to the 14ka-L simulation 
(Fig. 5e). The 14ka-5G results in a greater contribution of 
Antarctic bottom water into the deep Atlantic. In the Pacific, 
14ka-5G shows an increased contribution of Antarctic inter- 
mediate water, but the differences do not extend to include 
Antarctic bottom water (Fig. 5f). 

Proxy records suggest that 14 ka AMOC strength was sim- 
ilar to that of the LGM (McManus et al., 2004; Lynch- 
Stieglitz et al., 2007). However, other records suggest a more 
intermediate rate of overturning in shallower waters between 
the relatively elevated values of the LGM and the reduced 
values in the Holocene (Gherardi et al., 2009). In the Pacific, 
overturning at 14ka was in the middle of the transition from 
glacial to modern conditions (Herguera et al., 2010), which 
is similar to our results. 
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4 Discussion 

4.1 The effects of ice-sheet topography 

Uncertainty in the height of the LIS has a measurable impact 
on the simulated climate in G1SS ModelE2-R, particularly 
at 21 ka, where global mean SATs are significantly differ- 
ent between 21ka-L and 21ka-5G. This response to changes 
in LIS topography alone is due to a series of differences in 
the 21 ka climate, initiated by the differences in atmospheric 
circulation that arise from the enhanced topographic barrier 
in the LIS of 21ka-5G. In both 21 ka simulations, the polar 
jet is forced to the south of the LIS (Fig. Sle in the Supple- 
ment). However, the elevated topographic barrier in the 21ka- 
5G LIS drives the polar jet to be more zonal, whereas the 
21ka-L jet circulation has a greater meridional component, 
reflecting a shift in the downstream stationary waves in the 
Northern Hemisphere relative to 21ka-5G, as seen in 500 mb 
height differences (Fig. 2f). The primary impact of these sta- 
tionary wave differences is downstream colder temperatures 
across Siberia, Beringia, and the North Pacific, leading to ele- 
vated sea ice in the Sea of Okhotsk and Bering Sea, as well as 
enhanced snow cover across Siberia (Fig. SI in the Supple- 
ment and Fig. 4). Both of these impacts lead to an increase in 
ground albedo, thus reducing the total incoming shortwave 
budget in 21ka-L (Fig. 2g and h). This snow-cover-albedo 
positive feedback thus leads to further global cooling. 

In the 14 ka simulations, the difference in the magnitude of 
maximum LIS height is smaller between the two reconstruc- 
tions, but there is still downstream cooling over Siberia and 
Beringia in 14ka-L related to a similar difference in atmo- 
spheric circulation and stationary waves (Fig. S2e and fin the 
Supplement). However, differences in snow cover and sur- 
face albedo are no longer present across this region (Fig. 2g), 
and the global mean temperatures between 14ka-L and 14ka- 
5G are equivalent. This comparison with the 21 ka simula- 
tion might suggest a minimum cutoff in LIS maximum eleva- 
tion difference that induces strong enough changes in atmo- 
spheric circulation to cause significant differences in global 
mean surface temperature through a snow-albedo feedback. 
However, the 14 ka simulations are also forced with elevated 
boreal summer insolation and greenhouse gases, such that 
generally warmer globally averaged temperatures at 14ka 
(relative to 21 ka) may preclude the expansion of snow and 
sea ice in this region, eliminating the possibility of induc- 
ing this snow-albedo cooling feedback, even with the differ- 
ences in planetary wave strength induced by LIS topographic 
differences. 

Significant differences in ocean circulation also arise due 
changes in the LIS topography. All of our simulations show 
increased AMOC transport relative to the control, with south- 
ward displacement of the North Atlantic deepwater forma- 
tion (Figs. S3e and S5e in the Supplement). Glacial sur- 
face winds are enhanced over the North Atlantic (relative 
to 0 ka) in relation to the southward shift in the polar jet 


(Figs. Sid and S2d in the Supplement). Such windier con- 
ditions may provide the mechanism driving enhanced over- 
turning in our simulations (Wunsch, 2003; Montoya and 
Levermann, 2008). In addition, the 21ka-L simulation shows 
a slightly elevated freshwater balance (P — E, Fig. 2b) over 
the eastern North Atlantic, consistent with lower SSS, which 
may also contribute to the weaker overturning circulation rel- 
ative to 21ka-5G. 

In each of the 21 and 14 ka simulation pairings, the higher 
LIS (21ka-5G and 14ka-5G) results in a stronger AMOC 
(Figs. 4e and 5e). In either case, there is no significant change 
in river runoff between the simulations, but both 21ka-5G 
and 14ka-5G have stronger wind speeds over the North At- 
lantic relative to their 21ka-L and 14ka-5G counterparts, 
leading to enhanced wind-driven overturning. These differ- 
ences in AMOC strength are transferred to the deep Pa- 
cific with enhanced contribution of Antarctic bottom and 
intermediate waters associated with the stronger AMOC 
(Chikamoto et al., 2012). Despite these changes in circula- 
tion, the 21ka-L simulation has colder ocean temperatures 
(relative to 21ka-5G) in both the globally averaged ocean 
and the deep ocean, suggesting ocean temperatures are re- 
flecting the overall colder climate and reduced heat content 
of the system, as would be expected from the snow-albedo 
feedback mechanism driving colder SAT. 

4.2 Implications for ice-sheet stability 

The extent and volume of the LIS was determined by a num- 
ber of glaciological factors that impact ice-sheet stability. 
LIS surface mass balance was paced by boreal summer inso- 
lation (Hays et al., 1976; Imbrie and Imbrie, 1980). Such con- 
tinuous cycling of boreal summer insolation limits the size 
of the LIS through latitudinal shifts in the equilibrium line 
altitude moving in step with insolation (Oerlemans, 1991; 
Ruddiman, 2001), with variation in latitudinal extent of the 
LIS likely leading to ice volume hysteresis (Abe-Ouchi et 
al., 2013). Ice-sheet size is also limited by an ice-thickness- 
basal-melting negative feedback, whereby increases in ice 
thickness can lead to a decrease in the pressure melting point, 
decoupling the ice sheet from its bed and allowing for en- 
hanced basal sliding and ice height drawdown (Clarke, 1987; 
Payne, 1995; Marshall and Clark, 2002). These thermal in- 
stabilities may be exacerbated by subglacial till rheology 
conditions that might also limit ice-sheet height (Clark, 1994; 
Clark and Pollard, 1998; Licciardi et al., 1998). 

Our results may suggest an additional limit on LIS height. 
As shown in our LGM simulations, the larger LIS in 21ka- 
5G leads to a warmer global mean SAT due to atmospheric 
circulation changes and the snow-albedo feedback in north- 
eastern Asia. That an increase in LIS maximum height 
leads to global surface warming might suggest a climat- 
ically driven negative feedback on LIS surface mass bal- 
ance that limits ice-sheet height above a certain threshold. 
Satellite -based gravity field measurements of glacial isostatic 
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adjustment have suggested the need for revisions to ICE-5G 
(Peltier and Drummond, 2008; Peltier, 2009) to be included 
in the upcoming ICE-6G that lower maximum LIS eleva- 
tion (Peltier, 2010; otherwise see PMIP3 boundary condi- 
tions, http://pmip3.lsce.ipsl.fr/). These revisions may imply 
that LIS elevation in 1CE-5G is above the threshold of this 
elevation-warming feedback. In order to test this feedback 
mechanism, more analysis should be conducted on the im- 
pact of this SAT difference and ice-sheet topography on LIS 
surface mass balance. 

4.3 Implications for testing LGM climate 
reconstructions 

We present a range of climate responses to LGM external 
and internal forcings that develop solely due to two differ- 
ent reconstructions of the LIS. This range of uncertainty in 
this boundary condition alone is enough to provide signif- 
icant differences in LGM climate. Previous sensitivity ex- 
periments testing the impact of “ice” vs. “no ice” in model 
boundary conditions have shown large differences in LGM 
climate due to changes in atmospheric circulation and its at- 
tendant influence on ocean circulation (Justino et al., 2004; 
Otto-Bliesner et al., 2006; Abe Ouchi et al., 2007; Pausata et 
al., 2011; Hofer et al., 2012; Tharammal et al., 2012). Our 
21 ka model results show that even a 20 % change in LIS ele- 
vation between two LIS reconstructions has a similar impact 
on atmospheric and oceanic circulation and their influence 
on the LGM climate state. This suggests that the range of un- 
certainty within the existing LIS reconstructions can lead to 
significantly different results in simulated regional and global 
climate. 

The ability of GCMs to capture LGM AMOC is often 
used as a test of model skill (e.g., Timm and Timmerman, 
2007; Liu et al., 2009). Despite a range of proxy estimates 
for an AMOC target value (Yu et al., 1996; Adkins et al., 
2002; McManus et al., 2004; Curry and Oppo, 2005; Lynch- 
Steiglitz et al., 2007), coupled climate models simulate dif- 
ferent strengths and depths of the AMOC (Otto-Bliesner et 
al., 2007; Weber et al., 2007). The strength of the PMIP ap- 
proach toward simulating the LGM is in the assessment of 
climate variability (i.e., AMOC) across a variety of model 
physics and design with a common fixed set of boundary con- 
ditions. However, our results suggest that AMOC strength 
in a single model may vary by nearly 10% in response to 
changes in LIS topography alone, meaning that any uncer- 
tainties in the LIS boundary condition may be translated 
into the uncertainty in this important test of model skill. 
Thus, the uncertainty in the range of PMIP assessed AMOC 
strength might be expanded with inclusion of LIS topo- 
graphic uncertainty. 


4.4 Implications for climate sensitivity 

The LGM provides the most recent large-scale change in 
global climate and greenhouse gas concentrations through 
which global climate sensitivity can be assessed (Cruci- 
fix, 2006; Hansen et al., 2008; Schmittner et al., 2011; 
Hargreaves et al., 2012; PALAEOSENS, 2012). However, 
our model results suggest that significant differences in 
global mean temperature arise due to variability in LIS to- 
pography alone, indicating that the use of models to con- 
strain CO 2 sensitivity for the LGM should include some as- 
sessment of this boundary condition uncertainty in the over- 
all range of possible sensitivity estimates. This analysis is 
particularly relevant in the assessment of how uncertainty in 
boundary conditions leads to the development of fast feed- 
backs (i.e., snow-albedo) that are important in driving sen- 
sitivity (PALAEOSENS, 2012). Our 14 ka simulations show 
that even small differences in LIS height can lead to differ- 
ences in atmospheric circulation, but these differences are not 
enough to initiate the snow-albedo feedback. Depending on 
the model (or other boundary conditions), this precondition- 
ing may be sufficient to initiate fast-feedback mechanisms 
that lead to significant differences in global mean SAT, while 
other models may not provide the same response. Changes in 
LIS height within the uncertainty of reconstruction estimates 
may have a similar effect in other models. 

Given the large uncertainties in the high-latitude boundary 
conditions and their associated impact on high-latitude cli- 
mate, Hargreaves et al. (2012) instead correlated simulated 
LGM SAT from 20° S to 30° N with the global mean SAT 
change from CO 2 doubling sensitivity (2 x CO 2 ) simulations 
of each model in a grouping of PMIP2-CMIP3 pairings, 
which circumvented earlier issues with correlating global 
LGM SAT with 2 x CO 2 simulations (Crucifix, 2006). Our 
results suggest an additional source of high-latitude variabil- 
ity in LGM SAT that arises due to uncertainty in LIS topog- 
raphy alone, giving credence to the former study’s focus on 
tropical SAT at the LGM to constrain CO 2 sensitivity. How- 
ever, this additional source of uncertainty from ice-sheet to- 
pography in a region where LGM SAT anomalies are already 
at their greatest also suggests that focusing on the tropics 
could underestimate the full range of uncertainty in CO 2 sen- 
sitivity estimates. 

4.5 Uncertainty in other boundary conditions 

Finally, we have only assessed the impact of changing LIS 
height in one GCM. A number of other LGM boundary con- 
ditions have a significant impact on the global mean climate 
state, particularly dust/aerosol forcing (Penner et al., 2004; 
Lohmann and Feichter, 2005; Forster et al., 2007; Chylek 
and Lohmann, 2008) and vegetation land cover (Jahn et al., 
2005). The large range in uncertainties of the LGM and 
deglacial dust forcing (Petit et al., 1981; Thompson et al., 
1995; Mahowald et al., 1999; Reader et al., 1999; Mahowald 
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et al., 2006a, b; Power et al., 2008) and vegetation dynam- 
ics (Jackson et al., 2000; Prentice and Jolly, 2000; Harrison 
et al., 2001; Ray and Adams, 2001; Bigelow et al., 2003; 
Williams, 2003; Pickett et al., 2004; Williams et al., 2011) 
suggests the need to test the sensitivity of the climate re- 
sponse to the uncertainties in these boundary conditions. In 
the case of G1SS ModelE2-R, future sensitivity studies will 
have offline coupling of ModelE2 boundary conditions to the 
Lund-Potsdam-Jena (LPJ) dynamic vegetation model (Sitch 
et al., 2003; Gerten et al., 2004; Bondeau et al., 2007) to as- 
sess the impact of this land surface choice and any associated 
internal feedbacks in the model. 


5 Conclusions 

We attempt to assess the range of climate uncertainty that 
results from variability in the possible reconstructions of the 
LIS elevation at two different time periods, 21 and 14 ka. The 
simulated climate at each time slice results in significant dif- 
ferences in atmospheric and oceanic climate. The 21ka-L is 
significantly colder than the 21ka-5G simulation in both SAT 
and ocean temperatures, which is due to a snow-albedo feed- 
back in northeastern Asia that reduces the shortwave contri- 
bution to the system. This enhanced snow cover over north- 
eastern Asia in 21ka-L develops due to a diminished atmo- 
spheric ridge relative to 21ka-5G. The higher elevation LIS 
in 21ka-5G is a larger impediment to atmospheric circula- 
tion, impacting the location of the atmospheric jet through 
resulting differences in downstream planetary wave genera- 
tion, as has been suggested previously with LGM sensitivity 
tests removing the LIS (Justino et al., 2004; Otto-Bliesner 
et al., 2006; Abe Ouchi et al., 2007; Pausata et al., 2011; 
Hofer et al., 2012; Tharammal et al., 2012). We conclude 
that the range of uncertainty in existing LIS reconstructions 
is enough to have a similar “planetary wave effect” on North- 
ern Hemisphere atmospheric circulation, contributing to the 
large differences in temperature, snow cover, and surface 
albedo over northeastern Asia. These fully coupled AOGCM 
results build upon earlier lower-resolution simulations that 
indicate a similar dynamical adjustment with the comparison 
of ICE-4G and ICE-5G LIS reconstructions in an Earth sys- 
tem model of intermediate complexity (Justino et al., 2006), 
suggesting a robust impact of paleotopographic uncertainty 
on reconstructed climates in a variety of models. Future re- 
search should work on minimizing the uncertainty in the LIS 
reconstructions, thus reducing their impact on the uncertainty 
in simulated glacial climate. As we show in our 21 ka simula- 
tions, this LIS topographic uncertainty may significantly im- 
pact estimates of LGM climate sensitivity through the impact 
on global mean SAT. 

Given the significant differences between the resulting cli- 
mates from changing LIS topography alone, it is tempting to 
use climate proxy data to test which reconstruction is closer 
to the actual LGM LIS. However, the regions with largest 


model differences are the regions with the least proxy cov- 
erage and with the greatest variation between proxy recon- 
structions. In particular, the region of Siberia, Beringia, and 
the North Pacific has the largest modeled differences in SAT, 
<5 18 0, surface albedo, SST, and sea concentrations, all due to 
variation in LIS topography alone in both 21 and 14 ka sim- 
ulations. In theory, this region could serve as an important 
model-data test location on the LIS reconstructions used in 
GCMs. Unfortunately, the distribution of LGM SST proxy 
records from the North Pacific is limited (Kucera et al., 2005; 
Waelbroeck et al., 2009) and most records that do exist are 
confined along coastal regions, where GCMs may not ade- 
quately resolve changes in the Kuroshio and California Cur- 
rents, and sea level boundary conditions may also signifi- 
cantly influence ocean current behavior (Ortiz et al., 1997; 
Kao et al., 2006). Similarly, few terrestrial reconstructions of 
SAT and <$ 18 O a exist from lake records in the region (Melles 
et al., 2012; Peterse et al., 2011; Viau et al., 2008; Kurek et 
al., 2009; Bartlein et al., 2011), but perennial ice cover, proxy 
uncertainties, and the possibility of no-analog environments 
in pollen reconstructions may also limit the ability of such re- 
constructions to distinguish the climate differences discussed 
in this paper. To better constrain the topographic uncertainty 
of the LIS through the comparison of simulated climates in 
GCMs, more data are needed from regions where resulting 
climate differences are the greatest. Future data collection 
should focus on northeastern Asia and the North Pacific to 
help test LIS boundary conditions. 

In addition to large differences in atmospheric circula- 
tion and attendant impacts on surface conditions, the two 
LIS simulations at each time slice result in significant dif- 
ferences in the AMOC, with the larger ICE-5G ice sheet 
leading to stronger overturning at each time slice, suggest- 
ing that any given model’s ability to capture LGM anoma- 
lies in ocean overturning may be influenced as much by the 
LIS boundary condition as by limitations in model physics. 
At present, proxy uncertainty may preclude determination of 
whether the LGM AMOC was different than for the modern 
period (Burke et al., 2011). This uncertainty in LGM ocean 
circulation and surface temperature, particularly in regions 
where model differences are maximized, limits the determi- 
nation as to which LIS reconstruction is more correct in our 
simulations. 

One metric that may help differentiate between the LIS 
reconstructions is the Kr/ N 2 reconstruction of global mean 
ocean temperature (Headly and Severinghaus, 2007), as it 
should be indicative of the overall thermal state of the global 
system. Here, the 21ka-L simulation results in a global mean 
ocean temperature change more in line with the proxy re- 
constructions, whereas the 21ka-5G mean ocean is too warm 
(see Sect. 3.2.1). As this oceanic temperature difference is 
the result of the different atmospheric cooling from changes 
in planetary albedo, we woidd suggest that a lower eleva- 
tion LIS topography than 5G may be closer to what ex- 
isted at the LGM. However, the uncertainty in Kr/N 2 global 
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mean ocean temperature is large relative to the reconstructed 
change, therefore limiting the distinction of the two LGM 
simulations based on this proxy. Further refinement of LGM 
and 14ka climate reconstructions (particularly from Siberia 
and the North Pacific) will aid in the indirect determination 
of LIS topography, just as refinement of LIS topographic re- 
constructions will reduce inherent model uncertainty due to 
the current range in LIS topographic boundary conditions. 

Supplementary material related to this article is 
available online at http://www.clim-past.net/10/487/2014/ 
cp-10-487-2014-supplement.zip. 
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